% ----------------------------------------------------------------
% LaTeX Lecture
%  Conics
% **** -----------------------------------------------------------
\documentclass[a4paper, 10pt]{article}

\usepackage{amsmath}
\usepackage{fancyhdr}
\usepackage{layout}
\addtolength{\hoffset}{-1.0cm} \addtolength{\textwidth}{2cm}
\addtolength{\voffset}{-1.0cm} \addtolength{\textheight}{2cm}
\usepackage{appendix} %生成附录
\usepackage{array}
\usepackage{cite}
\usepackage{makecell}%能画表格粗线\Xhline{1.5pt}
% \usepackage[T1]{fontenc}
% \usepackage{textcomp}
\makeatletter
\def\@cite#1#2{\textsuperscript{[{#1\if@tempswa , #2\fi}]}}
\makeatother

\usepackage{appendix} %生成附录
\usepackage{flushend, cuted} %
\usepackage[dvipdfmx, %,pdftex这里决定运行文件的方式不同
                  pdfstartview=FitH,
                  CJKbookmarks=true,
                  bookmarksnumbered=true,
                  bookmarksopen=true,
                  colorlinks, %注释掉此项则交叉引用为彩色边框(将colorlinks和pdfborder同时注释掉)
                  pdfborder=001, %注释掉此项则交叉引用为彩色边框
                  linkcolor=black,
                  anchorcolor=green,
                  citecolor=green
                  ]{hyperref}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \href{url}{text}                           %
% \url{url}                                  %
% \nolinkurl{url}                            %
% \hyperbaseurl{url}                         %
% \hyperimage{imageurl}{text}                %
% \hyperdef{category}{name}{text}            %
% \hyperref{url}{category}{name}{text}       %
% \hyperref[label]{text}                     %
% \hyperlink{name}{text}                     %
% \hypertarget{name}{text}                   %
% \phantomsection                            %
% \addcontentsline{toc}{chapter}{\indexname} %
% \printindex                                %
% \autoref{label}                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\usepackage{listings}

%\usepackage{indentfirst,latexsym,bm}
\usepackage{amsmath,amssymb,amsfonts}
\usepackage{pifont}
\usepackage{CJK}
\usepackage[rgb]{xcolor}
\usepackage[dvipdfmx]{graphicx}
\usepackage{subfigure}
\usepackage{url}
\usepackage{mathrsfs}
\usepackage{caption}
\usepackage{shapepar}
\usepackage{float}
\usepackage[all,pdf]{xy}%状态转移图
\usepackage{tikz}%复杂的图
\usepackage{siunitx}
\usepackage{minted}
\usepackage{mdframed}
\usepackage{multirow}
%\surroundwithmdframed{minted}
\setlength\mathsurround{3pt}
\setlength{\columnsep}{3em}          %设置分栏间隔
\setlength{\parindent}{2em}          %设置段首缩进量
\renewcommand{\baselinestretch}{1.3} %重设行距       \usepackage{graphicx}
\newcommand{\lishu}[1]{   \CJKfamily{li}   {#1}}   %隶书\renewcommand*\fuigurename{\KAI{图}}
\newcommand{\heiti}[1]{  \CJKfamily{hei}  {#1}}   %黑体
\newcommand{\kaiti}[1]{  \CJKfamily{kai}  {#1}}   %楷体\usepackage{cite}
%\newcommand{\red}[1]{  \textcolor{red}  {#1}}   %红色\makeatletter
%\newcommand{\blue}[1]{ \textcolor{blue} {#1}}   %蓝色\def\@cite#1#2{\textsuperscript{[{#1\if@tempswa , #2\fi}]}}
%\newcommand{\green}[1]{\textcolor{green}{#1}}   %绿色\makeatother

\newcommand{\KAI}[1]{\begin{CJK*}{UTF8}{kai}{#1}\end{CJK*}}
\newcommand{\HEI}[1]{\begin{CJK*}{UTF8}{hei}{#1}\end{CJK*}}
\newcommand{\LI}[1]{\begin{CJK*}{UTF8}{li}{#1}\end{CJK*}}
\newcommand{\SONG}[1]{\begin{CJK*}{UTF8}{song}{#1}\end{CJK*}}
\newcommand{\YOU}[1]{\begin{CJK*}{UTF8}{you}{#1}\end{CJK*}}
\newcommand{\FSONG}[1]{\begin{CJK*}{UTF8}{fsong}{#1}\end{CJK*}}
% ----------------------------------------------------------------
\vfuzz2pt % Don't report over-full v-boxes if over-edge is small
\hfuzz2pt % Don't report over-full h-boxes if over-edge is small


%\usepackage{geometry}
%\geometry{left=2cm,right=2cm}
% MATH -----------------------------------------------------------
\DeclareMathOperator{\diag}{diag}
\DeclareMathOperator{\rank}{rank}
\DeclareMathOperator{\vecm}{vec}
\DeclareMathOperator{\vecs}{vecs}

\newcommand{\mfloor}[1]{ \left\lfloor {#1} \right\rfloor }
\newcommand{\mpair}[2]{ \left\langle {#1}, {#2} \right\rangle}


\renewcommand{\bf}[1]{\mathbf{#1}}
%\renewcommand{\vec}[1]{\bm{#1}}    %向量,  黑斜体
\newcommand{\mat}[1]{\bm{#1}}    %矩阵
\newcommand{\dif}{\mathrm{d}}
\newcommand{\diver}{\mathrm{div}}
\newcommand{\me} {\mathrm{e}}
\newcommand{\mi} {\mathrm{i}}
\renewcommand{\vec}[1]{\boldsymbol{#1}}
\newcommand{\vecmat}[1]{\vecm{\left( #1 \right)}}
\newcommand{\vecsmat}[1]{\vecs{\left( #1 \right)}}
\newcommand{\vecasym}[1]{[#1]_\times}   % antisymmetric matrix from a vector
\newcommand{\id} {\mathbbm{1}}   % identity operator
\newcommand{\fracode}[2]{\frac{\dif {#1}}{\dif {#2}}}         % ordinary differential operator
\newcommand{\fracpde}[2]{\frac{\partial {#1}}{\partial {#2}}} % partial differential operator
\newcommand{\fracpderow}[2]{\partial {#1}/\partial {#2}}
\newcommand{\fracoderow}[2]{\dif {#1}/\dif {#2}}
\newcommand{\fracpdemix}[3]{\frac{\partial^2 {#1}}{\partial {#2} \partial {#3}}}
\newcommand{\lap}[2]{\frac{\partial^2 {#1}}{\partial {#2}^2}}
\newcommand{\laprow}[2]{\partial^2 {#1}/\partial {#2}^2}
\newcommand{\secode}[2]{\frac{\dif^2 {#1}}{\dif {#2}^2}}
\newcommand{\set}[1]{\left\{ #1 \right\}}
\newcommand{\abs}[1]{\left| #1 \right|}
\newcommand{\absvec}[1]{\left| \bf{#1} \right|}
\newcommand{\ket}[1]{|#1 \rangle}
\newcommand{\bra}[1]{\langle #1 |}
\newcommand{\braket}[2]{ \langle #1 | #2 \rangle}
\newcommand{\norm}[1]{\lVert #1 \rVert}
\newcommand{\normF}[1]{{\parallel #1 \parallel}_\textrm{F}}
\newcommand{\tA}[1]{{#1}^\textsf{T}}
\newcommand{\inv}[1]{#1^{-1}}
\newcommand{\ginv}[1]{#1^+}    % Moore-Penrose (general) inverse
\newcommand{\tinv}[1]{{#1}^{-\textsf{T}}}
\newcommand{\linux}{$\mathrm{Linux}$}
\newcommand{\T}{^\textsf{T}}


\newcommand{\ES}[3]{\mathbb{#1}^{{#2}\times {#3}}}               % Euclidean space
\newcommand{\PS}[3]{\mathbb{#1}^{{#2}\times{#3}}}      % projective space
\newcommand{\dt}{\mathrm{d}t}
\newcommand{\hxq}{h_x^2}
\newcommand{\hyq}{h_y^2}
%----------------------------------------------------------------

\newcommand{\en}[1]{$\mathrm{#1}$}
\begin{document}
\begin{CJK*}{UTF8}{song}



%%---------------------------------
%使用 Scale=0.5的图片而且点击图片标题能够回到引用图片的地方, 需要和下面的\IREF命令配合使用(似乎和\par 有冲突, 多图时别在前面用\par )
\newcommand{\IMGD}[3]{
\begin{figure}[H]
\centering\includegraphics[scale=0.5]{#1}
\hyperref[#3totext]{\caption{#2}\label{#3}}
\end{figure}
}
%%---------------------------------
%IMaGe
\newcommand{\IMG}[4]{
\begin{figure}[H]
\centering\includegraphics[scale=#1]{#2}
\hyperref[#4totext]{\caption{#3}\label{#4}}
\end{figure}
}
%%---------------------------------
%使用点击文字到图片的超链接
\newcommand{\IREF}[1]{\hyperref[#1]{\KAI{图\ref{#1}}} \label{#1totext}}
%使用 Scale=0.5的图片而且点击图片标题能够回到引用图片的地方, 需要和下面的\IREF命令配合使用(似乎和\par 有冲突, 多图时别在前面用\par )

\newcommand{\TREF}[1]{\hyperref[#1]{\KAI{表\ref{#1}}} \label{#1totext}}


\newcommand{\FCZ}[1]{\begin{equation}#1\end{equation}}

%----------------------------------
%使用点击脚注到文字的超链接
\newcommand{\myfoot}[2]{\label{#1}\footnote{\hyperref[#1]{#2}}}

%%----------------------------------
%让小标题点击后能回到目录
\newcommand{\BC}[1]{\hyperref[contents]{#1}}

%%---------------------------------
%%同时设置工作报告的编号和标题
%\newcommand{\GZBG}[1]{\fancyhead[RO,LE]{WR-ABSR-03-No.#1-WZH}\title{\heiti{Weekly Research Report\\ WR-ABSR-03-No.#1}}}

%%---------------------------------
%%带圈的文字
\newcommand{\quan}[1]{\Large\textcircled{\small{#1}}}
%%---------------------------------
%% 设置页眉与页脚
%%*********************************************
\newtheorem{theorem}{定理}[section]
\pagestyle{fancy}
\fancyhead[RO,LE]{ES21-2015-GrA-P2011-64-WZH} % clear all fields
   %  请设置正确的个人文档编号
%\fancyhead[LO,RE]{\CJKfamily{hei>} \bfseries \leftmark}
\fancyhead[LO,RE]{一种\en{Laplace}方程的数值解法及程序实现}
%\fancyfoot[LE,RO]{\thepage}
\fancyfoot[RO,LE]{Academy of Blue Sky}
\fancyfoot[LO,RE]{SIAE}
\renewcommand{\headrulewidth}{0.4pt}
\renewcommand{\footrulewidth}{0.4pt}


% THEOREMS -------------------------------------------------------
\newtheorem{thm}{\textcolor[rgb]{0.98,0.00,0.00}{定理}}
\newtheorem{cor}[thm]{\textcolor[rgb]{0.00,0.00,1.00}{推论}}
\newtheorem{lem}[thm]{\textcolor[rgb]{0.00,0.00,1.00}{引理}}
\newtheorem{prop}[thm]{\textcolor[rgb]{0.00,0.00,1.00}{性质}}
\newtheorem{prob}[thm]{\textcolor[rgb]{0.98,0.00,0.00}{问题}}
\newtheorem{conj}[thm]{\textcolor[rgb]{1.00,0.00,1.00}{猜想}}
\newtheorem{defn}[thm]{\textcolor[rgb]{0.44,0.00,0.87}{定义}}
\newtheorem{rem}[thm]{\textcolor[rgb]{1.00,0.50,0.25}{注记}}
\newtheorem{eg}[thm]{\textcolor[rgb]{0.00,0.00,1.00}{例}}%[chapter]
\renewcommand{\abstractname}{摘要}
\renewcommand{\contentsname}{目录}
\renewcommand{\figurename}{\kaiti{图}}
%\renewcommand{\tablename}{\kaiti{表}}
\renewcommand{\refname}{\kaiti{参考文献}}
%\renewcommand{\refname}{\en{Reference}}

\setcounter{figure}{0}

%\fancyhead[RO,LE]{WR-ABSR-03-No.044-WZH}
\title{\heiti{一种\en{Laplace}方程的数值解法及程序实现}}


%\GZBG{081}
\author{王子昊\\110441328}

\date{\today}


%%*************************************************
%%  打印 标题, 作者, 日期等内容
%%*************************************************
\pagenumbering{roman}
\maketitle%

\begin{abstract}
本文提出了一种通用的\en{Laplace}方程的数值解法, 并使用Python作为编程语言将这种数值解法用程序实现. 在后期, 对这种算法做了进一步的优化, 使前者和后者的运行时间之比在$O(n^2)$量级. 最后, 探究了算法实现的时间复杂度问题, 并使用实验数据将我的优化算法的实现与SOR算法的实现法进行时间复杂度的比较.
\par
\vspace{1em}
\noindent\textbf{关键字：} \en{Laplace}方程; 数值解法; Python; 时间复杂度
\end{abstract}
\newpage
\tableofcontents

\label{contents}

%%*************************************************
%% 显示内容目录
%%*************************************************
\newpage
\pagenumbering{arabic}
%%*************************************************
%% 正文部分
%%*************************************************

\section{问题描述}
本文提出的数值解法适用于, 拥有第一类边界条件的, 矩形区域的, Laplace方程. 数学描述如下:

 $$ \displaystyle\frac{\partial^{2} u}{\partial x^{2}} +\frac{\partial^{2} u}{\partial y^{2}} = 0,\quad u|_{\partial\Omega}=f(x, y)$$
  $$\Omega: \left\{ \begin{array}{c}
  x\in[x_{min}, x_{max}]\\ y\in[y_{min}, y_{max}]
\end{array}   \right.$$

\section{理论分析}
\subsection{差分格式}
在本文提出的解法中, 使用的是传统的五点差分格式近似, 其数学描述如下:
在$x$方向上,$y$方向上均使用二阶中心差:
\begin{align}
\label{secondcenterx}
  (\frac{\partial^{2} u}{\partial x^{2}})_{i,j} &= \frac{\partial^{2} u_{i,j}}{\partial x^{2}} \approx \frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{h_x^2} \\
\label{secondcentery}
  (\frac{\partial^{2} u}{\partial y^{2}})_{i,j} &= \frac{\partial^{2} u_{i,j}}{\partial y^{2}} \approx \frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{h_y^2}
\end{align}
  $$h_x = (x_{max}-x_{min})/(N_x+1),\quad h_y = (y_{max}-y_{min})/(N_y+1)$$
其中, $N_x$, $N_y$分别为$x$方向上,$y$方向上的新增的节点数. 将式\eqref{secondcenterx},\eqref{secondcentery}代入\en{Laplace}方程, 有:
\begin{align}
\label{fivepoint}
\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{h_x^2} + \frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{h_y^2} = 0
\end{align}

\subsection{迭代方程的推导}
\subsubsection{代数方程矩阵化}
为了推导方便, 将等式\eqref{fivepoint}两边同时乘以${h_x^2}{h_y^2}$, 合并同类项后, 式\eqref{fivepoint}变为:
\begin{align}
	\label{fivepointex}
\hyq u_{i-1,j}+\hyq u_{i+1,j} + \hxq u_{i,j-1}+\hxq u_{i,j+1}-2(\hxq+\hyq)u_{i,j} = 0
\end{align}
对于每一个节点$(i,j)$都可以列出一个这样的方程, 将已知的边界节点项移至等式右边, 将所有节点的方程都这样列出即可得到这样的一个方程组. 按列选取节点构造向量$\vec{u}$:
$$\vec{u}=[u_{_{1,1}},u_{_{1,2}}\cdots u_{_{1,N_y}},u_{_{2,1}}\cdots u_{_{N_x,1}} \cdots u_{_{N_x,N_y}}]\T$$
通过仔细观察, 可将方程组转化为矩阵方程的形式:
\begin{equation}\label{matrixeq}
\vec{A}\vec{u}=\vec{b}
\end{equation}
其中,$\dim(\vec{A})=(N_x*N_y, N_x*N_y)$, $\vec{A}$有如下形式:
\begin{equation}\label{matrixa}
	A =\left[
		\begin{array}{ccccccc}
			\vec{D} & \vec{H} & \vec{Z} & \cdots &\cdots & \cdots & \vec{Z} \\
			\vec{H} & \vec{D} & \vec{H} & \ddots &\ddots & \ddots & \vdots \\
			\vec{Z} & \vec{H} & \vec{D} & \ddots &\ddots & \ddots & \vdots \\
			\vdots & \ddots& \ddots & \ddots & \ddots & \ddots & \vdots \\
			\vdots & \ddots& \ddots & \ddots & \vec{D} & \vec{H} & \vec{Z} \\
			\vdots & \ddots& \ddots  & \ddots & \vec{H} & \vec{D} & \vec{H} \\
			\vec{Z} & \cdots &\cdots & \cdots & \vec{Z} & \vec{H} & \vec{D} \\
		\end{array}
	\right]
\end{equation}
其中$\dim(\vec{D})=\dim(\vec{H})=\dim(\vec{Z})=(N_y, N_y)$, $\vec{D},\vec{H},\vec{Z}$的具体形式如下:
\begin{align*}
	D = \left(
		\begin{array}{cccccc}
			-2(\hxq+\hyq) & \hxq & 0 & \cdot & \cdot & 0 \\
			\hxq & -2(\hxq+\hyq) & \hxq & 0 & \cdot & 0 \\
			0 & \hxq & -2(\hxq+\hyq) & \hxq & 0 & \cdot \\
			0 & \cdot & 0 & \hxq & -2(\hxq+\hyq) & \hxq \\
			0 & \cdot & \cdot & 0 & \hxq & -2(\hxq+\hyq) \\
		\end{array}
	\right)
\end{align*}
\begin{align*}
	H = \left[
		\begin{array}{cccccc}
			\hyq & 0 & 0 & \cdot & \cdot & 0 \\
			0 & \hyq & 0 & 0 & \cdot & 0 \\
			0 & 0 & \hyq & 0 & 0 & \cdot \\
			0 & \cdot & 0 & 0 & 0 & 0 \\
			0 & \cdot & 0 & 0 & 0 & 0 \\
			0 & \cdot & \cdot & 0 & 0 & \hyq \\
		\end{array}
	\right], \quad
		Z = \left[
		\begin{array}{cccccc}
			0 & 0 & 0 & \cdot & \cdot & 0 \\
			0 & 0 & 0 & 0 & \cdot & 0 \\
			0 & 0 & 0 & 0 & 0 & \cdot \\
			0 & \cdot & 0 & 0 & 0 & 0 \\
			0 & \cdot & 0 & 0 & 0 & 0 \\
			0 & \cdot & \cdot & 0 & 0 & 0 \\
		\end{array}
	\right]
\end{align*}
\subsubsection{迭代矩阵}
将$\vec{A}$分解为
\begin{equation}
  \vec{A} = \vec{M} - \vec{N}
\end{equation}
带入式\eqref{matrixeq}中整理得
\begin{equation}\label{spliteq}
  \vec{u} = \vec{M}^{-1} \vec{N} \vec{u} + \vec{M}^{-1} \vec{b}
\end{equation}
式\eqref{spliteq}中$\vec{M}^{-1} \vec{N}$就是迭代矩阵. 将$\vec{A}$不同的分解诞生出了各种不同的方法, 原则上, 只要$\vec{M}$非奇异即可, 但事实上, 一般选取求逆容易的矩阵. 经过分析$\vec{A}$的结构以后, 将$\vec{M}$设为由$\vec{A}$的对角线组成的矩阵, 由于$\vec{A}$对角线上的元素都一样, 因此$$\vec{M}=A_{0,0}\vec{I}$$ 
从而有
$$\vec{M}^{-1}=\frac{1}{A_{0,0}}\vec{I},\quad \vec{N}=A_{0,0}\vec{I}-\vec{A}$$ 
将式\eqref{spliteq}离散化后即可得到迭代方程:
\begin{equation}\label{itereq}
  \vec{u}^{(k+1)} = (\vec{I}-\frac{1}{A_{0,0}}\vec{I}\vec{A}) \vec{u}^{(k)} + \frac{1}{A_{0,0}}\vec{I}\vec{b}
\end{equation}
\section{程序实现}
使用Python作为编程语言, 编程实现共分为以下五步:
\begin{itemize}
\item 确定边界条件
\item 生成式\eqref{matrixeq}中$\vec{A}$和$\vec{b}$
\item 利用$\vec{A}$和$\vec{b}$生成式\eqref{itereq}中的迭代矩阵
\item 根据式\eqref{itereq}迭代循环
\item 输出结果并画图
\end{itemize}
迭代循环的步骤如下:
\begin{itemize}
\item $\vec{u}^{(k)}$通过式\eqref{itereq}得到$\vec{u}^{(k+1)}$
\item 对$\vec{u}^{(k)}$和$\vec{u}^{(k+1)}$做差的结果取无穷范数, 得到误差
\item 判断误差是否满足精度要求, 满足:返回结果, 不满足:使用$\vec{u}^{(k+1)}$继续迭代
\end{itemize}
\section{运行结果}
\begin{figure}[H]
\includegraphics[scale=5.8]{RESULT.eps}
\includegraphics[scale=5.8]{ERROR.eps}
\caption{给定参数的运行结果图及其误差}
\end{figure}
%\IMG{6}{RESULT.eps}{给定参数的运行结果图}{result}
%\IMG{6}{ERROR.eps}{上图的误差}{error}
\section{优化}
\subsection{性能分析}
在能够输出正确结果后, 我发现, 运算时间对于$\vec{A}$的维数的变化非常敏感, 当$\dim(\vec{D})=(30,30)$的时候, 计算时间为1.09306秒, 而当$\dim(\vec{D})=(60,60)$时候, 计算时间上涨到了50.366538秒. 对此我很是不满意, 因此对Python程序进行了性能分析, 结果如下:
\begin{verbatim}
Total time: 51.3007 s
File: Laplace_Solver.py
Timer unit: 1e-06 s
Function: Method_Iterative at line 9

Line#  Hits     Time   Per Hit  %Time  Line Contents
===============================================================================
   9                                  @profile
  10                                  def Method_Iterative(A, b, Precision):
  11                                     """Solve \vec{A}\vec{x}=\vec{b}"""
  12     1     118449 118449.0    0.2  	  Iteration_Matrix, Constant_Matrix,\ 
                                              Initial_Vector =\
  13                                          Generate_Matrix_for_Iterate(A, b)
  14     1          3      3.0    0.0  	  Error = Precision + 1
  15     1          1      1.0    0.0  	  X = Initial_Vector
  16     1          0      0.0    0.0  	  Iteration_Times = 0
  17  4084       3807      0.9    0.0  	  while Error > Precision:
  18  4083   50854934  12455.3   99.1  	    	  X_next = Iteration_Matrix*X\
                                             	     + Constant_Matrix
  19  4083     315903     77.4    0.6  	    	  Error = np.max(np.abs(X_next - X))
  20  4083       4742      1.2    0.0  	    	  X = X_next
  21  4083       2898      0.7    0.0  	    	  Iteration_Times += 1
  22     1          0      0.0    0.0  	  return X, Error, Iteration_Times
===============================================================================
\end{verbatim}
根据上表第18行的数据可以得出, 超过99\%的时间都用在了迭代中计算即第18行的计算式子. 简单分析就可以知道, 消耗最多时间的是$(N_x*N_y,N_x*N_y)$维的矩阵与$(N_x*N_y,1)$维的列向量相乘, 简单分析可知, 单步消耗时间会是$O(N_x^2*N_y^2)$, 同时迭代步数会随着$N_x$和$N_y$的增加线性增加, 为了简化计算不妨设$N=N_x=N_y$, 总体计算的时间复杂度为$O(N^5)$.
\par
对于这种算法时间复杂度的不满意使我开始研究如何能够让我的程序运行的更快, 最先想到的是减少迭代次数. 在不改变整体结构的前提下想要找到减少迭代次数的方法, 我尝试了艾特肯加速收敛方法. 但遗憾的是, 程序实现时我发现这种方法对并不收敛, 因此放弃了. 
\subsection{对于矩阵乘法的优化}
\subsubsection{思想}
在后来的一段时间当中我就在思考如何简化矩阵的乘法. 在无意中, 我发现, 程序实现时, 单位阵与普通矩阵相乘以及零矩阵与普通矩阵相乘都需要花去一定量的时间, 而且这个时间随着矩阵维数的增加增加的很快, 也就是说, 计算机在做这些乘法的时候, 没有对单位阵和零矩阵做优化. 
\par
知道这一点后, 我做的第一件事就是把我程序中所有有关单位阵的乘法都略去. 然后我就在思考, $\vec{A}$可以表示为分块矩阵, 而且这个矩阵中有很多很多$0$, 那么矩阵$\vec{I}-\frac{1}{A_{0,0}}\vec{I}\vec{A}$是不是也可以写成有很多$0$的分块矩阵呢? 如果可以的话, 我可以把$(\vec{I}-\frac{1}{A_{0,0}}\vec{I}\vec{A}) \vec{u}_{k}$拆分开来, 只将$\vec{I}-\frac{1}{A_{0,0}}\vec{I}\vec{A}$中非零分块矩阵与$\vec{u}_{k}$中相应元素相乘即可, 这样可以省去大量与零相乘的运算. 同时, 这样还有另外一个好处, 就是不必计算以及储存$\vec{A}$矩阵, 只需计算和储存组成$\vec{A}$矩阵的分块矩阵$\vec{D},\vec{H}$即可, 大大节省了储存空间. 
\subsubsection{数学表达}
将式\eqref{itereq}中的迭代矩阵$\vec{I}-\frac{1}{A_{0,0}}\vec{I}\vec{A}$分块化:
\begin{equation}\label{matrixdprime}
\vec{I}-\frac{1}{A_{0,0}}\vec{I}\vec{A}=
\left(
		\begin{array}{ccccccc}
			\vec{D'} & \vec{H'} & \vec{Z} & \cdots &\cdots & \cdots & \vec{Z} \\
			\vec{H'} & \vec{D'} & \vec{H'} & \ddots &\ddots & \ddots & \vdots \\
			\vec{Z} & \vec{H'} & \vec{D'} & \ddots &\ddots & \ddots & \vdots \\
			\vdots & \ddots& \ddots & \ddots & \ddots & \ddots & \vdots \\
			\vdots & \ddots& \ddots & \ddots & \vec{D'} & \vec{H'} & \vec{Z} \\
			\vdots & \ddots& \ddots  & \ddots & \vec{H'} & \vec{D'} & \vec{H'} \\
			\vec{Z} & \cdots &\cdots & \cdots & \vec{Z} & \vec{H'} & \vec{D'} \\
		\end{array}
	\right),\quad \begin{array}{c}
					\vec{D'}=\vec{I}-\frac{1}{A_{0,0}}\vec{D}\\
					\\
					\vec{H'}=-\frac{1}{A_{0,0}}\vec{H}
					\end{array}
\end{equation}
将列向量$\vec{u}$分块化:
$$\vec{u}^{(k)}=
\left[
	\begin{array}{c}
		u_{_{1,1}}^{(k)}\\
		u_{_{1,2}}^{(k)}\\
		\vdots\\
		u_{_{1,N_y}}^{(k)}\\
		u_{_{2,1}}^{(k)}\\
		\vdots\\
		u_{_{N_x,1}}^{(k)}\\
		\vdots\\
		u_{_{N_x,N_y}^{(k)}}
	\end{array}
\right] = 
\left[
	\begin{array}{c}
		\vec{U}_{_{1}}^{(k)}\\
		\vec{U}_{_{2}}^{(k)}\\
		\vdots\\
		\vdots\\
		\vdots\\
		\vec{U}_{_{N_x}}^{(k)}
	\end{array}
\right],\quad 
\vec{U}_{_{n}}^{(k)} = 
\left[
	\begin{array}{c}
		u_{_{n,1}}^{(k)}\\
		u_{_{n,2}}^{(k)}\\
		\vdots\\
		u_{_{n,N_y}}^{(k)}
	\end{array}
\right]
$$
从而, 若假设$N=N_x=N_y$, 则矩阵乘法由$(N*N,N*N)*(N*N,1)$变为$2*(N,2*N)*(2*N,1)+(N-2)*(N,3*N)*(3*N,1)$, 即
\begin{equation}\label{itermeq}
\centering
\begin{split}
\vec{U}_{_{1}}^{(k+1)}&=
\left[
		\begin{array}{cc}
			\vec{D'} & \vec{H'} 
		\end{array}
\right]
\left[
	\begin{array}{c}
		\vec{U}_{_{1}}^{(k)}\\
		\vec{U}_{_{2}}^{(k)}
	\end{array}
\right]\\
\vec{U}_{_{n}}^{(k+1)}=&
\left[
		\begin{array}{ccc}
			\vec{H'} & \vec{D'} & \vec{H'} 
		\end{array}
\right]
\left[
		\begin{array}{c}
			\vec{U}_{_{n-1}}^{(k)}\\
			\vec{U}_{_{n}}^{(k)}\\
			\vec{U}_{_{n+1}}^{(k)}
		\end{array}
\right]\\
\vec{U}_{_{N_x}}^{(k+1)}&=
\left[
		\begin{array}{cc}
			\vec{H'} & \vec{D'} 
		\end{array}
\right]
\left[
	\begin{array}{c}
		\vec{U}_{_{N_x-1}}^{(k)}\\
		\vec{U}_{_{N_x}}^{(k)}
	\end{array}
\right]
\end{split}
\end{equation}
其中, $n\in[2,N_x-1]$.
\subsubsection{程序实现}
同样使用Python作为编程语言, 编程实现共分为以下五步:
\begin{itemize}
\item 确定边界条件
\item 生成式\eqref{matrixa}中的$\vec{D}$和$\vec{H}$
\item 利用$\vec{D}$和$\vec{H}$生成式\eqref{matrixdprime}中的$\vec{D'}$和$\vec{H'}$
\item 根据式\eqref{itermeq}进行循环迭代
\item 输出结果并画图
\end{itemize}
迭代循环的步骤如下:
\begin{itemize}
\item 将$\vec{u}_{_{k}}$拆分成$\vec{U}^{_{k}}$依次代入式\eqref{itermeq}得到所有的$\vec{U}^{(k+1)}$
\item 将所有的$\vec{U}^{(k+1)}$拼成$\vec{u}^{_{k+1}}$
\item 对$\vec{u}^{(k)}$和$\vec{u}^{(k+1)}$做差的结果取无穷范数, 得到误差
\item 判断误差是否满足精度要求, 满足:返回结果, 不满足:使用$\vec{u}^{(k+1)}$继续迭代
\end{itemize}
\subsubsection{结果验证}
对于同样的边界条件, 分别统计使用原始算法和对矩阵乘法优化后算法计算$X$方向上节点数从$4-100$的情况下的计算时间, 这里为了达到最大的计算量, 设$Y$方向上节点数与$X$方向上节点数相同, 因此, 总结点数就是从$16-10000$, 纵轴时间单位为秒, 结果如\IREF{omt}所示.
\IMG{0.5}{TimeCost.eps}{原始算法和矩阵乘法优化后算法的比较}{omt}
\subsection{对于迭代的优化}
\subsubsection{思想}
在开始着手写这篇文档的时候我还没有想到这个方法, 但是写着写着我突然想起老师曾经说过的, 在一次迭代过程中, 将已经生成的一部分结果放入本次迭代过程的余下步骤中作为已知量进行迭代, 能够起到减少迭代次数的作用, 有时候还能够起到保持辛特性的作用. 于是放下文档的编写开始实现这种想法.
\subsubsection{数学描述}
分块化的迭代矩阵由式\eqref{itermeq}, 变为如下形式:
\begin{align}
\label{itermatrealtimeeq1}
\vec{U}_{_{1}}^{(k+1)}&=
\left[
		\begin{array}{cc}
			\vec{D'} & \vec{H'} 
		\end{array}
\right]
\left[
	\begin{array}{c}
		\vec{U}_{_{1}}^{(k)}\\
		\vec{U}_{_{2}}^{(k)}
	\end{array}
\right]\\
\label{itermatrealtimeeq2}
\vec{U}_{_{n}}^{(k+1)}=&
\left[
		\begin{array}{ccc}
			\vec{H'} & \vec{D'} & \vec{H'} 
		\end{array}
\right]
\left[
		\begin{array}{c}
			\vec{U}_{_{n-1}}^{(k+1)}\\
			\vec{U}_{_{n}}^{(k)}\\
			\vec{U}_{_{n+1}}^{(k)}
		\end{array}
\right]\\
\label{itermatrealtimeeq3}
\vec{U}_{_{N_x}}^{(k+1)}&=
\left[
		\begin{array}{cc}
			\vec{H'} & \vec{D'} 
		\end{array}
\right]
\left[
	\begin{array}{c}
		\vec{U}_{_{N_x-1}}^{(k+1)}\\
		\vec{U}_{_{N_x}}^{(k)}
	\end{array}
\right]
\end{align}
其中, $n\in[2,N_x-1]$. 与式\eqref{itermeq}的区别在于, 式\eqref{itermatrealtimeeq2}与式\eqref{itermatrealtimeeq3}中右侧分块向量中的第一个分量是由第$k+1$次迭代得到的.
\subsubsection{程序实现}
同样使用Python作为编程语言, 编程实现共分为五步, 与对乘法优化后的步骤相同, 不同的是循环迭代的步骤:
\begin{itemize}
\item 将$\vec{u}^{_{k}}$深度复制一份设为$\vec{v}^{{k}}$
\item 将$\vec{u}^{{k}}$拆分成$\vec{U}^{{k}}$代入式\eqref{itermeq}得到$\vec{U}^{(k+1)}$
\item 利用得到的$\vec{U}^{(k+1)}$更新$\vec{u}^{{k}}$中对应的那部分元素
\item 完成$\vec{u}^{{k}}$的全部更新, 将其记为$\vec{u}^{{k+1}}$
\item 对$\vec{u}^{(k+1)}$和$\vec{v}^{(k)}$做差的结果取无穷范数, 得到误差
\item 判断误差是否满足精度要求, 满足:返回结果, 不满足:使用$\vec{u}^{(k+1)}$继续迭代
\end{itemize}
\subsubsection{结果验证}
对于同样的边界条件和网格节点划分条件, 分别统计使用矩阵乘法优化算法和对矩阵乘法优化后再对迭代优化的算法的计算时间, 纵轴时间单位为秒, 结果如\IREF{oit}所示.
\IMG{0.5}{TimeCost1.eps}{矩阵乘法优化算法与矩阵乘法优化+迭代优化算法的时间比较}{oit}
\par
同时这两种算法的迭代次数比较如\IREF{oii}所示.
\IMG{0.5}{IterationTimes.eps}{矩阵乘法优化算法与矩阵乘法优化+迭代优化算法的迭代次数比较}{oii}
\section{时间复杂度的探究}
\subsection{思想}
如何衡量一种算法程序实现的时间复杂度? 我想是不是可以利用做实验拟合数据的方式得时间与算法中某一变量的关系? 由于之前做了很多实验, 因此手里有很多节点数与程序运行时间的数据, 因此我用这些数据来分析程序运行的时间与节点数的关系.
\subsection{程序实现与结果}
程序实现非常简单, 设$Y$方向上节点数与$X$方向上节点数相同即$N=N_x=N_y$, 将程序运行时间和单轴节点数$N$导入MATLAB, 调用cftool中的指数拟合($Time = a*Node^b$)进行拟合, 得到如下结果.
\IMG{0.6}{Original.eps}{原始算法的单轴节点数与运行时间的关系$T = 6.008^{-9}*N^{5.618}$}{originalnt}
\IMG{0.6}{Optimised.eps}{矩阵乘法优化算法的单轴节点数与运行时间的关系$TO = 1.425^{-6}*N^{3.757}$}{optmnt}
\par
简单分析可得优化算法节省的时间是接近$N^2$级别的($N^{1.861}$).
\par
最后简略的看下周露莎同学使用MATLAB实现的$SOR$方法\cite{bianweihui}的单轴节点数与运行时间的关系图.
\IMG{0.6}{SOR.eps}{$SOR$算法的单轴节点数与运行时间的关系$SOR_Time = 5.026^{-10}*Node^{6.538}$}{sor}
\par
可以看出, $SOR$算法的实现花费的时间比原始算法的实现花费的要多, 但这并不能说明$SOR$算法不够好, $SOR$算法能够显著的减少迭代步数, 但是对于寻找最优松弛因子的时候会有性能开销, 再加上MATLAB本身的在非计算环节的性能开销, 得到这个结果也是合理的. 但是究竟为什么会这样? 如何定量的分析一个算法实现的时间复杂度? 如何比较一个算法本身的时间复杂度与其实现的时间复杂度? 这些都是可以探讨的问题.
\section{结论}
本文提出了一种解Laplace方程的数值解法, 重点工作是对于这种方法进行了优化, 得到了不错的结果, 即在节点数为$100\times100$时时间节省10倍以上. 在结尾处对于如何衡量一种算法实现的时间复杂度提出了自己的方法.
\section{致谢}
感谢\KAI{张鸿燕}老师, 他给我提供了充足的学习资源, 用于解决这个问题所需的资料都在手边, 大大节省了时间. 同时感谢周露莎学姐, 我修改并运行了她写的SOR方法的MATLAB程序用于进行性能分析.

\addcontentsline{toc}{section}{参考文献}
\bibliographystyle{unsrt}
\bibliography{bibliography.bib}
\newpage
\addcontentsline{toc}{section}{附录: 程序代码}
%\section*{Appendix: Code}
\renewcommand\listoflistingscaption{附录: 程序代码}
%\renewcommand\listingscaption{Program code}
\usemintedstyle{tango}
\definecolor{bg}{rgb}{0.9,0.9,0.9}
\listoflistings
\begin{mdframed}[frametitle={File: Laplace\_Solver.py}, backgroundcolor=bg]
\begin{minted}[linenos=true, breaklines]{python3}
import numpy as np
from matplotlib.pyplot import plot, savefig, title, legend, ylim, cla, xlabel, ylabel, annotate, figure, show
from time import clock
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter

def Method_Iterative(A_shape, b, D, H, Precision):

	D_Iteration, H_Iteration, Constant_Matrix, Initial_Vector = Generate_Matrix_for_Iterate(A_shape, b, D, H)
	Ny = D_Iteration.shape[0]
	Iteration_Matrix_First = np.bmat([D_Iteration, H_Iteration])
	Iteration_Matrix_Middle = np.bmat([H_Iteration, D_Iteration, H_Iteration])
	Iteration_Matrix_Last = np.bmat([H_Iteration, D_Iteration])
	Nx = int(A_shape[0]/Ny)
	X_Real_Time = Initial_Vector
	Error = Precision + 1 
	count = 0
	while Error > Precision:
		X = np.copy(X_Real_Time)
		X_Real_Time[0:Ny] = Calcule_Part_of_X(Iteration_Matrix_First, Constant_Matrix[:Ny], X_Real_Time[:Ny*2])
		for i in range(1,Nx-1):
			X_Real_Time[Ny*i:Ny*(i+1)] = Calcule_Part_of_X(Iteration_Matrix_Middle, Constant_Matrix[Ny*i:Ny*(i+1)], X_Real_Time[Ny*(i-1):Ny*(i+2)])
		X_Real_Time[Ny*(Nx-1):] = Calcule_Part_of_X(Iteration_Matrix_Last, Constant_Matrix[Ny*(Nx-1):], X_Real_Time[Ny*(Nx-2):])
		# Error = np.linalg.norm(X_Real_Time - X)
		# Error = np.sum(np.abs(X_Real_Time - X))
		Error = np.max(np.abs(X_Real_Time - X))
		count = count + 1
	return X, Error, count

def Calcule_Part_of_X(Iteration_Matrix, Constant_Matrix, Initial_Vector):
	X = Iteration_Matrix*Initial_Vector + Constant_Matrix
	return X

def Generate_Matrix_for_Iterate(A_shape, b, D, H):
	"""
	A = M - N 
	N = M - A
	x = M^{-1}*N*x + M^{-1}*b
	x = Iteration_Matrix*x + Constant_Matrix
	"""
	"""
	把 Iteration_Matrix像对A一样分块， 这样可以分块的运算X
	"""
	Dimention = D.shape[0]
	M_inv = 1/D.item(0,0)
	D_Iteration = np.matrix(np.eye(Dimention,Dimention)) - M_inv*D
	H_Iteration = -M_inv*H
	Initial_Vector = np.matrix(np.zeros([A_shape[0],1]))
	Constant_Matrix = M_inv*b
	return D_Iteration, H_Iteration, Constant_Matrix, Initial_Vector


def Generate_Matrix_A_b(Nx, Ny, Xmax, Xmin, Ymax, Ymin):
	"""
	Nx和Ny分别表示了在X方向和Y方向上划分的点数，因而网格数需要加1
	"""
	dx = (Xmax - Xmin)/(Nx+1)
	dy = (Ymax - Ymin)/(Ny+1)
	D, H = Generate_Diagonal_Matrix(Nx, Ny, dx, dy)
	A_shape = (Nx*Ny, Nx*Ny)
	b = np.zeros(Nx*Ny)
	Border_Up, Border_Down, Border_Left, Border_Right = Border_Generator(Nx, Ny)
	for i in range(Nx):
		for j in range(Ny):
			if 0 == i:
				b[Ny*i+j] = b[Ny*i+j] + (-dy**2)*Border_Left[j+1]
			if (Nx-1) == i:
				b[Ny*i+j] = b[Ny*i+j] + (-dy**2)*Border_Right[j+1]
			if 0 == j:
				b[Ny*i+j] = b[Ny*i+j] + (-dx**2)*Border_Up[i+1]
			if (Ny-1) == j:
				b[Ny*i+j] = b[Ny*i+j] + (-dx**2)*Border_Down[i+1]
	b = np.mat(b).T
	return A_shape, b, D, H

def Generate_Diagonal_Matrix(Nx, Ny, dx, dy):
	Block_Line_String = str(dx**2)+","+str((-2)*(dx**2+dy**2))+","+str(dx**2)+",0"
	Block_Line_String = Block_Line_String + ",0"*(Ny-3)
	index = Block_Line_String.find(",")+1
	Diagonal_Matrix_String = Block_Line_String[index:]
	position = Block_Line_String.rfind(",")
	for i in range(Block_Line_String.count(",")-1):
		Diagonal_Matrix_String = Diagonal_Matrix_String+";"+"0,"*i +Block_Line_String[0:position]
		position = Block_Line_String.rfind(",",0,position)

	D = np.matrix(Diagonal_Matrix_String)
	Kernal_Dim = D.shape[0]
	H = np.matrix(np.eye(Kernal_Dim))*(dy**2)
	return D, H

def Border_Generator(Nx, Ny, Border_Type="array"):
	Border_Up = np.linspace(0,0,Nx+2)
	Border_Down = np.linspace(5,5,Nx+2)
	Border_Left = np.linspace(0,5,Ny+2)#+np.cos(np.linspace(0,5,Ny+2))
	Border_Right = np.linspace(0,5,Ny+2)#+np.sin(np.linspace(0,5,Ny+2))
	if Border_Type == "array":
		return Border_Up, Border_Down, Border_Left, Border_Right
	if Border_Type == "matrix":
		return np.mat(Border_Up[1:-1]), np.mat(Border_Down[1:-1]), np.mat(Border_Left).T, np.mat(Border_Right).T

def Plot_Result(Z, Nx, Ny, Xmax, Xmin, Ymax, Ymin):
	Border_Up, Border_Down, Border_Left, Border_Right = Border_Generator(Nx, Ny, Border_Type="matrix")
	xs = np.linspace(Xmin, Xmax, Nx+2)
	ys = np.linspace(Ymin, Ymax, Ny+2)
	X, Y = np.meshgrid(xs, ys)
	Z = np.reshape(Z, [Ny,Nx],"F")#reshape the matrix as Fortran-like
	Z = np.bmat([[Border_Left, np.bmat([[Border_Up], [Z], [Border_Down]]), Border_Right]])
	Z = np.array(Z)

	fig = figure()
	ax = fig.gca(projection='3d')

	cset = ax.contourf(X, Y, Z, zdir='z', offset=np.min(Z), cmap=cm.coolwarm)
	cset = ax.contourf(X, Y, Z, zdir='x', offset=-1, cmap=cm.coolwarm)
	cset = ax.contour(X, Y, Z, zdir='y', offset=Ymax*1.05, cmap=cm.coolwarm)
	surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, alpha=0.5, cmap=cm.coolwarm, linewidth=0.01)
	fig.colorbar(surf, shrink=0.5, aspect=5)
	#ax.zaxis.set_major_locator(LinearLocator(10))
	#ax.zaxis.set_major_formatter(FormatStrFormatter('%.1e'))
	ax.set_xlabel('X')
	ax.set_xlim(Xmin, Xmax*1.05)
	ax.set_ylabel('Y')
	ax.set_ylim(Ymin, Ymax*1.05)
	ax.set_zlabel('Z')
	ax.set_zlim(np.min(Z), np.max(Z)*1.05)
	show()

if __name__ == "__main__":

	Nx = 20
	Ny = 20
	A_shape, b, D, H = Generate_Matrix_A_b(Nx, Ny, 5, 0, 5, 0)
	start_time = clock()
	Z, Error, count = Method_Iterative(A_shape, b, D, H, 1e-5)
	end_time = clock()
	print("\nIterative Method:", "\n\tMatrix Size:", A_shape, "\n\tError:", Error, "\n\tIterative times:", count)
	print("\nUse time:", end_time - start_time, "s.")
	Plot_Result(Z, Nx, Ny, Xmax=5, Xmin=0, Ymax=5, Ymin=0)
\end{minted}
\end{mdframed}
\newpage
\end{CJK*}
\newpage
\end{document}
